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ABSTRACT 

-^ , We introduce an adaptation of the well known Tree+SPH numerical scheme to Post 

^^ I Newtonian (PN) hydrodynaniics and gravity. Our code solves the (0+1+2. 5)PN equa- 

^^ I tions. These equations include Newtonian hydrodynaniics and gravity (OPN), the first 

o ' order relativistic corrections to those (IPN) and the lowest order gravitational radia- 

•<^ ■ tion terms (2.5PN). We test various aspects of our code using analytically solvable test 

if^ • problems. We then proceed to study the IPN effects on binary neutron star coales- 
cence by comparing calculations with and without the IPN terms. We find that the 

^N^ , effect of the IPN terms is rather smah. The largest effect arises with a stiff equation of 

^ ', state for which the maximum rest mass density increases. This could induce black hole 

l/^ [ formation. The gravitational wave luminosity is also affected. 

Subject headings: gravitation,hydrodynamic,relativity,stars: neutron 

Q^' 1. Introduction 

O 

5-H ' The development of numerical methods for the solutions of full 3D general relativistic hydro- 

^ ■ dynamic problems is still at a preliminary stage (Font et al. 1998). In the meanwhile, various 

approximate approaches to this problem have been attempted. Of these some include an approx- 
imation to the metric in the form of the conformal flatness condition (CFC) (Mathews & Wilson 
1997) or by using a Post-Newtonian (PN) formulation of the equations (Shibata et al. 1997; Oohara 

5^ I &; Nakamura 1997). Recently it has been suggested that at least in some cases these two approxi- 

mations are of the same order of accuracy (Kley & Schafer 1999). In this work we set out to modify 
the popular Tree-SPH (Benz et al. 1990; Hernquist & Katz 1989) numerical scheme formulated 
for Newtonian self gravitating hydrodynamic systems to work in the PN approximation. We need 
to adapt both parts of the Tree-SPH scheme - the tree code and the SPH. Our goal is to study 
the coalescence of binary neutron stars (BNS), an intrinsically 3D hydro+gravity problem, and to 
begin the exploration of the general relativistic effects which are important in this process. 
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Smooth Particle Hydrodynamics (SPH) is a Lagrangian, grid-less particle niethod for solving 
the hydrodynamic equations. It is this lack of a grid which makes SPH especially appealing for 
the efficient solution of complex 3D problems. SPH has been used in hydrodynamic simulations 
including gravity, magnetic fields and in special relativistic problems (Kheyfets et al. 1990; Siegler 
&: Riffert 1999). Laguna et al. (1993) have also implemented a fully general relativistic SPH with 
a fixed Kerr background. For a review on SPH techniques see Monaghan (1992) and Benz (1990). 
The Barnes-Hut Tree algorithm (Barnes & Hut 1986) is a 0(iV log A^) method for calculating 
gravitational forces between N particles. It has been combined with SPH in order to get a powerful 
and efficient particle based gravity+hydro solver. 

We combine the efficiency of the SPH approach with a PN formalism suggested by Blanchet 
et al. (1990) (BDS). The PN approximation to gravity involves expanding the relativistic equations 
with a small parameter 0{v'^/c'^) where f is a typical velocity in the system (thermal velocity, 
orbital velocity etc...). IPN means we neglect all terms proportional to v* /c^ and higher. The IPN 
approximation includes in the first order, many relativistic effects. The main shortcoming of the 
IPN approximation is that it misses all effects which have to do with gravitational radiation. These 
effects appear only at the 2.5PN {y^ /e') order. The BDS formalism includes the Newtonian physics 
(OPN), the leading order PN effects (IPN) and the leading order gravitational radiation effects 
(2.5PN) in a self consistent way. It recasts the PN equations in a form similar to the Newtonian 
equations thus facilitating the adaptation of the Tree-SPH algorithm to solve this problem. 

In this work we introduce the code and present code tests. We also use the code to simulate 
binary neutron star (BNS) coalescence in the (0+1+2. 5)PN approximation and to compare these 
results to Newtonian simulations. In section 2 we introduce the numerical method. In section 3 we 
examine the results of various code tests. We present the results of the BNS coalescence simulation 
in section 4 and we conclude in section 5. 



2. Numerical method 

The BDS formalism recasts the (0+1+2. 5)PN gravitation and hydrodynamic equations in a 
form resembling the Newtonian (OPN) equations. This enables the solution of these equations 
using methods adapted from Newtonian gravity (such as the Tree+SPH method we use here). The 
formalism reduces all the relativistic non-local equations to compact supported Poisson equations. 
The PN order of the various terms in the equations of motion can be read off their coefficients 
- OPN terms have no coefficients, IPN have coefficients proportional to l/c? and 2.5PN terms 
have coefficients proportional to l/e' . This enables us to "turn off" various powers of the PN 
approximation by setting to zero the corresponding coefficients. We use this option later when 
making a (0+2.5)PN calculation. 

The independent matter variables used are the following set: p* the coordinate rest mass 
density, e* the coordinate specific internal energy and 'w the specific linear momentum. In fully 



relativistic terms these are defined as: 

P* = \/gu°P, (1) 

e* = e(p*), (2) 

Wi = {c^ + e+p/p) —, (3) 

where p is the rest mass density, e{p) the specific energy, p{e, p) the pressure and u^ the four- 
velocity (Greek indices run from O to 4, Latin indices from 1 to 3). The corresponding BDS 
variables are the above quantities neglecting ah terms except OPN, IPN and 2.5PN. Using these 
variables the formahsm yields an evolution system which consists of 9 Poisson equations and 8 
hyperbohc equations (see Appendix A). 

In the OPN approximation (c ^ oo) Wi = v^ , p^, is the mass density, the only needed potential 
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is U^:, and we recover the known Newtonian equations of motion. From Eq. (A7) we have Q- = 
d^Qij / dt^ + 0(1 / c^) , where Qij is the quadropole moment, enabhng us to compute the gravitational 
radiation luminosity of the system to lowest order as 

Low^feg'gg'. (4) 

The BDS equations are self-consistent as long as the system is only mildly relativistic. This 
can be characterized by noting that at least the parameters a/c^, (5/c^ and S/c? (Eq. Al, A2 
and A3 respectively) must be small. This means that contrary to OPN systems which are always 
self-consistent but can be physically wrong the IPN system will cease to be self consistent at the 
time when its results are too far from the general relativistic results. This inconsistency can be 
understood as follows: given two functions expanded in a small parameter rj <^\ 

f = fo + fiV + ---, (5) 

9 = go+9iV + ■■■, (6) 

the product of these functions is 

f 9 = fo9o + (/offi + fi9o)v + {fo92 + 2/151 + f29oW + ■■■■ (7) 

If we choose to truncate the functions at the first order (setting /« = O for i > 1) and indeed 77 ^ 1 
then we can also neglect the r/^ terms in the product, but if r/ ~ 1, we will also need to take into 
consideration the 'ifigirf term in the product if we wish to be self-consistent (the approximation 
would break down in any case giving wrong results). Similarly in order to make the IPN system self 
consistent we would need to add some 2PN terms to the equations complicating them even more. 
This is not the case if we wish to truncate the functions at the zero order (/j = O for i > 0). Then 
we can consistently truncate the product also at the zero order. This explains the self-consistency 
of the Newtonian (OPN) approximation. 



For the solution of the Poisson equations we used a Barnes-Hut tree (Barnes & Hut 1986). 
The Barnes-Hut tree is a 0{N log N) method commonly used for A^-body gravitational force cal- 
culations. The principle behind the Barnes-Hut tree is that at a given position, the gravitational 
force from a cluster of distant particles can be approxiniated by using only global properties of the 
cluster - the monopole, dipole and quadrupole moments of the niass distribution. The Barnes-Hut 
tree provides a simple way to calculate these moments and to cluster the particles. It is remarkable 
that the same formahsm can be used with minor modification to solve any compact supported 
elhptic equation. We adapt this method to solve a general compact supported ehiptic equation 
(see Appendix B) and we use it here to solve numerically equations (A8)-(A12). In all our runs 
we set the tree accuracy parameter 6 to 0.5. Gravitational softening was naturally incorporated by 
assigning particles with their respective SPH mass density in the gravitational calculations. 

The evolution equations were solved using SPH, a Lagrangian particle based scheme for solving 
hydrodynamic problems. The use of SPH was facilitated by the similarity of the equations in the 
BDS formalism to the Newtonian equations (indeed this is the whole point in the BDS formalism). 
The particle mass used in Newtonian SPH was replaced by the conserved mass m* = j d?xp^: in 
our scheme. When using SPH Eq. A17 is automatically satisfied and only equations A18-A20 have 
to be solved. The use of SPH requires adding some artificial viscosity in order to resolve shocks. We 
use the standard artificial viscosity (e. g. Monaghan 1992; Benz 1990) consisting of a term analogous 
to bulk viscosity and a Von Neuman-Richtmyer artificial viscosity term. 

We used a symmetrical form of the SPH equations which guarantees the exact conservation 
of the "momentum" J d?x p^wi in the hydrodynamics section of the code. The momentum isn't 
exactly conserved in the gravity section of the code because the tree algorithm introduces a small 
asymmetry into the forces. This problem arises also in Newtonian SPH calculations as well and it 
leads to spurious accelerations of the center of mass of the system. These acceleration are small for a 
typical Tree parameters (of the order of 10~^ compared to other accelerations in the system) and can 
be corrected by simply subtracting the center of mass acceleration from the particle accelerations 
at each time step. 

For time integration we used a second order Runge-Kutta integrator with adaptive step size 
control. The step size was determined so that none of the integrated variables will change by more 
than a predetermined amount, set to 0.5% in all runs. In addition we check that the Courant 
condition is always satisfied. We did not implement single particle time steps. The simulations 
were run on standard Pentium II 450MHZ workstations and took about 255 hours of CPU time 
(~ 11 days) each. 



3. Code Tests 

Numerous studies have been conducted on the capabilities and limitations of the SPH formalism 
(Hernquist & Katz 1989; Davies et al. 1993; Lombardi et al. 1999). We have put our code through a 



suite of tests to insure that our implementation is correct and that it reproduces previously obtained 
results using available codes. Specifically we have compared the Newtonian (OPN) version to the 
Newtonian SPH code used in Davies et al. (1994). This comparison ensures the validity of the 
Newtonian part of our code. We devised other simple tests for the different PN aspects of the code 
- (0+l)PN hydrodynaniics, gravitational radiation damping (2.5PN) and (0+l)PN hydrodynamics 
and gravity. These all have an analytical or easily obtainable result and take a reasonable time to 
run. 

(0+l)PN hydrodynamics (without gravity) were tested using the ID shock tube problem. We 
have compared the results to both the Newtonian and the relativistic solutions (Taub 1948; McKee 
& Colgate 1973; Hawley et al. 1984). These tests were conducted with a ID version of our code 
which differs from the 3D version used elsewhere in this paper only by the normalization of the 
SPH smoothing kernel. The results of the tests are presented in Figure 1. These tests demonstrate 
that our PN code gives better results than a Newtonian code for mildly relativistic problems, 
with shock velocities of ~ 0.2c. These conditions are at the upper limit of the validity of the 
(0+l)PN approximation as e ~ O.lc^. We compare (0+l)PN calculation to the analytic Newtonian 
solution and the relativistic analytical solution. The error in the (0+l)PN velocity is smaller by 
an order of magnitude than the Newtonian error. This leads to a more accurate estimate of the 
shock's position. The errors for the energy density and mass density are larger, but still better or 
comparable to the Newtonian error. We conclude that even in these extreme conditions for the 
(0+l)PN hydrodynamic approximation, it fares at least as well as Newtonian hydrodynamics, and 
is much better at estimating the velocity of the fluid. In all the quantities, the relative error of 
the (0+1) PN result as compared to the relativistic analytical solution is of the order of 1% except 
for single particles which reside at discontinuities. In Figure (1) we also show the convergence of 
the results. We show the relative error compared to the analytic relativistic results for different 
particle numbers. Except at the discontinuities the error relative to the (0+l)PN result converges 
to ~ 10~^. The size of the zone around the discontinuities with a large error decreases when 
the number of particles increases, indicating, again, that the discontinuities affect only a fixed 
small (~ 10), number of particles. The convergence of the error to ~ 10~^ is consistent with the 
order of magnitude of the largest 2PN term - e^/c*^ which is the expected error of the (0+l)PN 
approximation. We point out that 3D tests of shocks with the particles positioned randomly, which 
are expected to model the conditions in a 3D run more realistically, give much poorer resolution 
(Rasio &: Shapiro 1991a). 

Gravitational radiation backreaction (the 2.5PN terms) was tested using two point masses in 
orbit. In this run, all the IPN terms were discarded leaving only (0+2.5)PN terms. Since we used 
two point masses, hydrodynamics didn't play any role in this test. The rate of change of the orbit's 
radius was compared to the analytical result of d/a oc — a~^ (see e. g. Shapiro & Teukolsky (1983), 
Chapter 16). Figure 2 depicts the results of this test. The two point masses were set initially in a 
Keplerian orbit. Since this is only an approximation to actual quasi-stationary orbits of (0+2.5)PN 
gravitation, a brief period of relaxation preceded the orbital decay after which the expected power 



law enierged. 

The (0+l)PN gravity and hydrodynamics were tested using a static polytropic star with a 
polytropic index 7 = 5/3. Initial conditions were constructed using a (0+l)PN expansion of the 
Oppenheinier-Volkoff (OV) equations for a relativistic spherical star. We ran this test with two 
resolutions, 2151 and 5140 particles respectively, for approximately 50 hydrodynamic time scales 
(which we take to be y/BF/GM). for the initial conditions we use regularly placed unequal mass 
particles. The stars were constructed by placing the particles in a face-centred cubic (FCC) lattice, 
thus niaxiniizing the number of neighbors each particle has. Using an iterative procedure, each 
particle is than assigned a mass so that the resulting density profile matches the OV density 
profile. 

Figure 3a depicts the radius enclosing 95% of the mass. This radius oscillates with a decreasing 
amplitude converging to a constant value. There are less than 50 oscillation in the graphs since the 
oscillation period for these stars is some factor of order unity times \J B? jGM^ which we take to be 
the hydrodynamical timescale. The oscillations are caused by several factors in the SPH algorithm. 
Every SPH particle adjusts it's size h so that it will have approxiniately 50 other particles as 
neighbors (within a sphere of radius 2/i). In our initial conditions all the particles have the same 
/i, this causes the particles near the star's surface to have only half the neighbors, and their h 
is increased by the code. This in turn lowers the density and takes the star out of equilibrium. 
Another reason for the oscillations is that VP 7^ O on the boundary of the star due to numerical 
errors. Finally in the initial conditions the particles are positioned on a Cartesian grid which clearly 
conflicts with the spherical symmetry of the OV density profile assigned to them. The star begins 
to oscillate but the oscillations are damped by the artificial viscosity present in the SPH algorithm. 
The expected errors in all quantities (radius, central density, etc.) from the SPH discretisation 
error {{h) /Ro)"^ are 1.5% and 0.9% for the 2151 and 5140 particle runs respectively. As can be seen 
in Figure 3, the oscillations decrease when we increase the resolution. The final rest mass density, 
shown in Figure 3b, converges towards the analytical curve as the resolution increases. 



4. PN BNS coalescence 

BNS's provide an excellent test-bed for gravitational and nuclear astrophysics. A binary pair 
of neutron stars will loose angular momentum and energy via gravitational wave emission as has 
been observed (see Taylor 1994, and references therein). This process will ultimately lead to a 
coalescence of the two neutron stars. The gravitational waves emitted from the coalescence are 
expected to be observed by gravitational wave detectors coming on-line in the next decade, such 
as LIGO (Barish 1998), VIRGO (Fidecaro & Collaboration 1997), GEO (Hough et al. 1996) and 
TAMA (Kozai k Tama-300 CoUaboration 1999). 

At the final stages of the coalescence the BNS system must be described using detailed 3D 
modeling of gravitational and hydrodynamic effects. This restricts the study of the last stages 



of coalescence to numerical niethods. Many groups have performed numerical simulations using 
difFerent approxiniations to this problem. Results have been obtained using Newtonian dynamics 
by Davies et al. (1994), Rasio & Shapiro (1995), Wang & Swesty (1997), Ruffert & Janka (1998) 
and Rosswog et al. (1999). Post Newtonian (PN) results have also been obtained by Shibata et al. 
(1997) and Oohara & Nakamura (1997). Recently an almost fully relativistic result was obtained 
by Mathews & Wilson (1997) who used the CFC approximation (although see Flanagan (1999) and 
Kley & Schafer (1999) for cautionary remarks on the validity of this approximation) . 

Although sophisticated techniques were developed in order to obtain equilibriuni binary con- 
figurations both in the (0+l)PN (Shibata 1997, 1998) and general relativistic (Baumgarte et al. 
1998) cases, converting the resulting density field into SPH particles is not a trivial task. Previous 
works using Newtonian SPH usually manufacture equilibrium configurations by relaxing the sys- 
tem in the co-rotating frame, in which the stars are stationary (e. g. Rasio & Shapiro 1991b). This 
approach was unavailable to us because of the complicated nature of the (0+1+2. 5)PN equations. 
Instead we assign two spherical stars with Keplerian velocity (as in Shibata et al. 1998), and after 
some oscillations the system settles down into a stationary state. 

We model the binary NS system by two equal polytropes with zero spins relative to an inertial 
observer (For further discussion on the initial spins and their implications see Rosswog et al. (1999) 
and references therein). On top of this we made several other simplifying assumptions in our 
calculation which are already discussed in Davies et al. (1994), namely ignoring neutrino transport. 
The masses we use for each star are less than M0 for radii of about 30 Km. Although these 
parameters are far from realistic for NS's, they allow us to investigate the effect of general relativity 
on the coalescence while still in the regime of applicability of the (0+1+2. 5)PN approximation for 
the whole duration of the run. We compare the (0+1+2. 5)PN results (hereafter denoted P) to 
(0+2.5)PN runs (hereafter denoted N) with the same mass and initial separation to highlight the 
IPN effects. We make a total of 6 runs (3 pairs of P-N runs). The parameters for these different 
runs are summarized in Table 1. 



Run 


7 


M [Mo] 


R^ [Km] 


Particles 


mcsc [Mo] (# part.) 


P 


N 


P 


N 


1 


5/3 


0.52 


29.8 


33.9 


4996 


< 3.7 X 10-5 (0) 


< 3.7 X 10-5 (0) 


2 


2.6 


0.50 


34.7 


36.9 


20562 


2.6 X 10~^ (8) 


6.9 X 10-^ (20) 


3 


2.6 


0.92 


30.5 


35.6 


18520 


< 4.3 X 10-5 (0) 


4.5 X 10-4 (6) 



Table 1: Summary of results and parameters for the various runs. M is the conserved mass per star, 
Rif the radius encompassing 95% of the conserved mass of each star in the initial configuration, 
rriesc is the mass of unbound particles in the final configuration. The upper limits on nicsc are the 
masses of the lightest particles in the run. 



When comparing a P run to a N run with the same mass, we must take into account that a 
polytrope that is static in (0+l)PN gravitation is not static in OPN, Newtonian, gravitation. For 



the sanie mass, the Newtonian polytrope has a larger radius. This causes the initial polytrope to 
grow at the beginning of the N runs as demonstrated in Figure 4. Note that unhke Figure 3a 
where we show the radii of isolated starts, here we show the radii of the stars in the first timesteps 
of the coalescence. In order to differentiate between IPN effects and the effects of these differences 
in initial conditions, we try to present all results scaled by the appropriate mass and initial radius. 
For instance, the dynaniical orbital instability sets in at a/i?* ~ 3 (Rasio &: Shapiro 1995) where 
a is the binary separation. For a N run this is at a larger separation than for a P run with similar 
mass. This means that the P runs will have to evolve longer with gravitational radiation damping 
as the only mechanism for decreasing the separation. This behavior can be seen in Figure 5a. 

The scaling we use during the rest of the discussion is the following: distance is measured in 
units of i?*, luminosity is measured in units of C^ / e' {M / R^,)^ , rest mass density (/9^=) is measured 
in units of M/R^ and energy in units of the rest mass energy of the star Mc^. Time is measured in 
hydrodynamic time scales (R^/GM)^'"^ and is shifted so that for each run the minimal separation 
occurs at t = 0. i?* is the initial radius of the stars. M = J d^zp^, is the conserved mass for each 
run. We emphasize that in the P runs the conserved mass is not the same as the rest mass. 

We define the center of mass of each star to be the center of conserved mass of the particles 
belonging to the star at the initial time step. Since we use a particle based scheme we can follow 
particles throughout the simulation and use this definition even after the stars have touched. In 
Figure 5a we show the separation of the centers of mass of the two stars as a function of time. As 
can be seen, the runs start in a slightly elliptical orbit which slowly decays. Two distinct processes 
cause this decay. The first is the conversion of orbital angular momentum into stellar spins Jg via 
gravitational torques. This process is also present in Newtonian gravity. The second process is 
the emission of gravitational radiation, carrying with it angular momentum Jg^. This is a 2.5PN 
process having no Newtonian parallel. The dynaniical orbit instability sets in for both the P and N 
runs as the separation reaches Si?*, {t ~ —25). This causes a rapid plunge, and a merger in about 
one orbital period. In Figure 6 we show the ratio Jgw/Js up to the dynamical instability. In the 
runs N3 and P3 this ratio is close to unity near t = — 25. In all other runs it is apparent that the 
spin-up of the stars plays a more important role than gravitational radiation in causing the merger. 
The stars touch at t ~ —15 (when a/R^ ~ 2). The minimal separation at t = O is achieved when 
the cores, which hold most of the mass, have merged. At t > O the centers of mass "bounce" . The 
bounce is more pronounced for the P2, N2 and N3 runs because of the stiffer EoS used in these 
runs. In the P3 run however, IPN gravity is strong enough to counteract this stiffness and the 
bounce is comparable to that of the Pl and Nl runs. 

In Figure 5b we depict the gravitational radiation luminosity (Eq. 4) of the merger. The 
characteristics of the luminosity peak are similar in all the runs once we allow for the different 
initial radii of the stars. The P3 run however exhibits a pronounced second peak at t ~ 5 at about 
the same luminosity of the system at the last orbits before the merger. This second peak is absent 
in all other runs. This distinct feature is a result of the strong IPN gravity and stiff EoS of the P3 
run. 



In order to explain this distinct feature of the P3 run we show in Figure 7 the maximum rest 
mass density /?* for all the runs. The Pl and Nl runs exhibit a larger relative density because of 
their softer EoS. In all the runs we see a dip in the maximuni density at the time of the rapid infall 
caused by the orbital instability. This corresponds to the stars shedding each other's mass as they 
move closer together. This stage ends at t « —15, when the stars touch, and is followed by a fast 
rise up to t ~ —7. The difference between the N2 and N3 runs and the Nl run can be attributed 
to the softer EoS of the latter while the difference between them and their respective P runs is due 
to the stronger gravitational attraction in the P runs. The maximum density in the N2 and N3 
runs does not have the peak at i ~ —7 which is evident in all other runs, most distinctly in the 
P3 run where it rises to about 10% more than it's final value. It is in this run where the stiff EoS 
and the strong IPN gravity combine to induce a large compression of the cores which delays their 
final merger into a single axisymmetric central object. This delay turns the merger into a two part 
process and produces the second peak in the luminosity at t ~ 3. This can be seen in Figure 8 
where we compare the cores of the P3 and N3 runs at t ~ —10,2,20. t ~ —10, is just after the 
first peak in the luminosity. We see the in the N3 run the cores have almost merged. At t « 2 the 
cores in the N3 run have merged and formed an almost axisymmetric object which emits very little 
gravitational radiation, on the other hand the P3 cores are still almost separate and certainly far 
from axial symmetry. At t ~ 20 the cores of the P3 run have already merged completely and emit 
little gravitational radiation as well. 

In this context it is interesting to note Ruffert et al. (1997) where there is a comparison of the 
gravitational wave luminosity produced by different numerical schemes. SPH is found to inhibit 
the second peak in the gravitational wave luminosity possibly because of the numerical viscosity. 
Using other numerical schemes the second peak in the luminosity is about one half the height of the 
first peak. This raises the possibility that IPN terms, may in reality have an even more prominent 
effect on the luminosity. 

The actual gravitational waveforms emitted by the systems are shown in Figure 9. Here 
we see a difference in the period of the waveforms between the P and N runs corresponding to a 
different orbital period before the actual merger. During and after the merger there is no qualitative 
difference between the waveforms. 

In Figure 10 we compare the total energy emitted by gravitational waves. We start the com- 
parison at t = —25 when all runs have roughly the same relative separation, before the stars touch. 
When comparing the similar mass runs Pl, Nl, P2, and N2 we see that a softer EoS implies more 
energy emitted in gravitational waves, while the more massive runs P3 and N3 emit an order of 
magnitude more energy as compared to run P2 and N2 which have a similar EoS. The P3 run emits 
almost twice the energy of the N3 run. We note that in this figure the axis scaling is different from 
Figure 5b e. g. run P3 and N3 have almost similar L in Figure 5b, but since it is scaled by R^ the 
actual P3 luminosity is a factor of ~ 2 higher than the N3 luminosity 

We now turn to look at the morphological differences between the runs. In Figures 11 and 
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12 we show the contours of coordinate rest mass density on the orbital (x — y) plane at various 
times. The difference due to the different EoSs is the most striking. The Pl and Nl runs lead to a 
final configuration with almost non-existent spiral arms, which are very prominent in the P2, N2, 
P3, and N3 runs. There is almost no difference between the Nl and Pl results, while we see that 
N2 and N3 result in longer spiral arms than P2 and P3 respectively. This is most prominent in 
the difference between N3 and P3. The central object at t = 20 is approximately the same size in 
all the runs and is axisyninietric. For a closer look at this central object we show contours of the 
rest mass density at t=20 in the x — z plane (Fig. 13). In all the runs we see a central core, with 
a density above O.IM/Rf and which has an equatorial radius of approximately 1.5i?*. The polar 
radius is smaller in the Nl and Pl runs than in the P2, N2, P3, and N3 runs. Surrounding this 
core is a halo with a radius of about 5i?*. The halo is extended vertically in the Nl and Pl runs 
up to a distance of 3i?*, while in the P2, N2, N3, and P3 runs it has a height of only 2i?^<. The 
P3 run resulted in the most compact object, as could be expected. In all runs, there is a funnel, a 
zone of low density, around the axis of rotation. 

Finally, we calculate how much mass can escape the system. We do this by counting all 
the particles in the final configuration which are at a distance greater than 6i?* from the origin 
and which have positive total Newtonian energy. The total Newtonian energy of a particle is 
-E'Newt = ^rnv'^ + me — mU^, with e the specific internal energy and U^, (Eq. A8) the Newtonian 
gravitational potential. Particles which are closer than QR^, will surely interact with other particles 
before exiting the system and so their energy won't be conserved. For particles further than this, 
the gravitational potential U^ is of the order of 0.01 and the velocity is of the order of O.lc which 
ensures that the Newtonian energy will be conserved. Only in run N2, P2 and N3 have any particles 
escaped, as shown in Table 1. The N2 has the highest escaped mass followed by N3. For the other 
runs we are only able to give upper bounds on the escaped mass by using the mass of the lightest 
particles in the run. It's these minimal mass particles which escape since they are located at the 
surface of the stars in the initial conditions. We note that escaped mass is also given in Rosswog 
et al. (1999). However a comparison between our results is impossible since we use different initial 
conditions and EoSs. 



5. Summary 

We have introduced here a PN adaptation of the Tree+SPH algorithm. This adaptation is 
made possible by the BDS formalism which recasts the (0+1+2. 5)PN equations of gravity and 
hydrodynamics in a form resembling the Newtonian equations. We have tested various aspects 
of our code by comparing the numerical results with known analytical solutions of relativistic 
problems. The (0+l)PN hydrodynamic part had been tested using a relativistic shock tube; the 
2.5PN gravitational radiation reaction terms were tested using two point masses in orbit and the 
IPN gravitation+hydrodynamic terms were tested using a spherical static OV polytrope. The code 
passes all these tests with the expected accuracy. 



11 



Using this code we have investigated the IPN effects on BNS coalescence. We compare runs 
with identical initial conditions but different physics. In the N runs, we use only the (0+2.5)PN 
terms, in the P runs we include also the IPN terms. In both runs we keep the 2.5PN gravitational 
radiation terms. These terms lead to a slow decreases of the orbital separation until a critical 
separations is reached. At this point a dynamical instabihty sets in and the stars merge within one 
orbital period. We use polytropes with a mass of less than Mq and a radius of about 30Km for 
the runs. These are not typical NS parameters but the (0+1+2. 5)PN approximation in the BDS 
formalism is valid only when the IPN terms are small compared to the OPN terms and this set an 
upper limit to the compactness of the stars. Therefore, Our results do not describe typical BNS 
coalescence, but rather the effect of the IPN terms on this process. 

Our results show that when going up to higher masses and thus to more relativistic conditions, 
there appears a prominent peak in the maximum rest mass density just before the cores merge in 
the P3 run. This peak is absent in the N3 run. This peak in rest mass density could mean a that 
the probability of the coalesced object collapsing into a black hole is larger than that estimated 
by Newtonian codes. Also, we see that the energy emitted in gravitational waves is almost twice 
as large in the P3 run as compared with the N3 run. This difference is also seen in the profile of 
the gravitational wave luminosity of the system. This result supports the suggested idea to use the 
gravitational wave signal as a probe on the details of the merger process. Our simulations show 
that the absence of a prominent second peak in the luminosity indicates a soft EoS. 



A. BDS Equations 

We use the Lagrangian formulation of the BDS equations (Blanchet et al. 1990). First we 
define some auxiliary quantities 

a = 2U^- 7, (iw^ + 3U,) , (Al) 

/3 = iw2 + 3t/* + £* + —, (A2) 

P* 

S = 3^2^e, +3^-[/,, (A3) 

p* 

Ai = AUi + \C,-\x'dsUs, (A4) 

^5 = l{R-Qf}x'djU,y (A5) 

Pij = 2 / cfix p^ i SwidjU^: — 2—djp^ + x'^WsdsjU^ — x^dsjUs 1 , (A6) 

where p* = p(e*,/9*) is the pressure, 7* = 9 log p* /9 log p* is the adiabatic index and w^ = WiWi. 
Using these quantities we solve the following Poisson equations 

A^J* = -AnGp^, (A8) 



12 



AC/, = 


- —AnGp^Wi, 


ACi = 


-- -AnGz'ds {p*Ws) , 


AU2 = 


= -4ttGp^6, 


AR = 


-- -47TGQfjx'djp,. 



(A9) 
(AlO) 

(Ali) 

(A12) 



The forces and the velocity are defined next 



1 1 4G ra 

v' = Wi-^Pwi + ^Ai + -^WsQfJ, (A13) 



c^ c^ 5 c^ 



^press ^ _]_g 



t 



1 + ^)^* 



(A14) 

p* L\ C"/ J 

i^'™ = (l + -^']d^U, + ^diU2-^Wsd^As, (A15) 

We are now in position to advance the system in time using the evolution equations 

p* = -p*div\ (A17) 

e; = -^div\ (A18) 

p* 

X' = v\ (A19) 

Wi = Ff''''' + Fl^^ + Fr^, (A20) 

where the dot represents the Lagrangian time derivative d = dta + v^diU. 



B. Adapting the Barnes-Hut tree to solve general Poisson equations 

The Barnes-Hut tree is a method for calculating the gravitational force between N particles 
in A^logA^ operations. The force on a given particle is calculated by summing the forces from 
individual particles if they are close, and by using a niultipole approximation for far away clusters 
of particles. A cluster is considered far away if d/l < 9 where d is the clusters size, / is it's distance 
and 9 is an external parameter governing the degree of approximation of the method. The tree 
itself is a data structure specially adapted for efficient clustering of particles and calculations of the 
multipole moments. In practice the multipole approximation is accurate enough to be stopped at 
the quadrupole term. 

Although it was originally invented for the calculation of gravitational forces, the Barnes-Hut 
tree can solve any compact supported Poisson equation with little modification. Given the multipole 
moments of the source of the equation, the potential and its derivatives can be calculated in an 
analogous way to the calculation of the gravitational potential and force in the following way: 
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Given a compact source rj and the masses and densities of each particle, we will solve the 
Poisson equation 

AU = -Attt]. (BI) 

The first three niultipole nioments of a cluster of particles are 

M^ = Y.—^v^ (B2) 

Di = E— ^p4' (B3) 

P Pp 

Qij = 2_^ ^P K^p^p ~ S^ij-'^p-^p) ' (^4j 

p Pp 

where the index p runs over the cluster particles, nip and pp are the particle mass and density 
(so that nip/pp ~ Vp the particle volume) rjp is the value of r] on the particle and Xp is the i-th 
component of the particle coordinate. 

Using these nioments the cluster's potential and it's derivatives can be calculated as follows: 
U" = Ml^£!^ + ^^ (B5) 

1 D'^Xi 

^j^jk 5XiXjQ^- 

+ 5 ~ o 7 ^fc 

dlkU'^ = :r-6lk + 3—^XkXl 

1 D'^Xi 

+3— {Dlxi + Dj'zfe - Dlxi5ik) + 15— y^Xfca;i 

where r is the distance to the cluster, 5ij is the Kronecker delta and a double index is summed 
over. 
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Fig. 1. — Shock tube results for a 7 = 5/3 gas. (a) The velocity. (b) The internal energy density 
e. (c) The rest mass density p. Initial conditions were p^i = 1, p„- = 0.05, and e* = 0.05 using 
800 particles. On the left, each graph shows the Newtonian analytical results in a solid line, the 
relativistic analytical results in a dashed line and the (O-l-l)PN numerical points as crosses. Under 
each graph is the error in the (O-l-l)PN result as compared to the relativistic result. On the right 
the each graph shows the error relative to the analytic relativistic solution. The vertical lines mark 
the left edge of the rarefraction, the right edge of the rarefraction, the contact discontinuity and 
the shock going from left to right. The error is calculated for solutions with 400, 800 and 1600 
particles. Except around the discontinuities the error converges to ~ 10~^, which is the expected 
error resulting from neglecting 2PN and higher ternis (e^ ~ TO^^). 
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Fig. 2. — Results of (0+2.5)PN runs with point masses which check the accuracy of gravitational 
radiation damping: a plot of —a/ a vs. a. The initial points with a > 0.7ao (where ao is the initial 
separation) were excluded to ahow for numerical relaxation. The expected slope is -4. 



18 



1.05 



0.99. 




O 



10 



20 3 1^1^ 

t [(R^/GM)^^^] 



40 



50 



xlO 



13 




Fig. 3. — (0+l)PN 7 = 5/3 polytropic stars made out of 2151 and 5140 particles. The simulation 
ran for about 50 hydrodynamic time scales. The parameters of the star were Rq = 34 Km and 
M = 0.49 M0. (a) The radius inclosing 95% of the mass. The sohd hne and dashcd hne are for 
the 2151 and 5140 particle runs respectively. (b) Rest mass density profiles. The initial profile in a 
sohd hne and the final profiles of the 2151 and 5140 particle runs in crosses and circles respectively. 
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Fig. 4. — The radii of the stars as a function of time for the beginning of the Pl and Nl runs. The 
initial polytropes are static in the (0+l)PN approxiination. The initial expansion in the N run 
reflects in the larger size of the N static polytrope. Based on this graph we take the radius of the 
Pl run to be 29.8Km and the radius of the Nl run to be 33.9Km. For the Pl run this is within 3% 
of the initial radius and for the Nl run this is 10% more than the initial radius. Note that in this 
graph the units are different than in all following graphs 
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Fig. 5. — (a) The separation between the centers of mass of the stars (in units of initial radius) as 
a function of time for all the runs. (b) The gravitational wave luminosity of the system for all the 
runs. 
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Fig. 6. — The ratio Jgw/Js of angular monientum in the gravitational radiation to angular momen- 
tum in stellar spins. The ration is shown up to the time of the dynamical orbital instability. In all 
runs but N3 and P3 spin-up of the stars is the dominant process for causing the merger. 
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Fig. 7. — Maxiinuiii coordinate rest mass density p* in the P and N runs. (a) the Pl and Nl runs. 
(b) The P2 and N2 runs. (c) The P3 and N3 runs. The maximum drops during the merger and 
chmbs back up afterwards. Notice the distinct rise at t ~ — 5 in the P3 run 
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Fig. 8. — Rest mass density p* of the cores of the N3 (top) and P3 (bottom) run. The contours are 
logarithmic with a spacing of lO*^'^. The length scale is in R^. The minimal contour is at 10% of 
the maximum rest mass density so the spiral arms are not visible in this plot. The cores is shown 
at times (a) t f» -10. (b) t ^2. (c) t ss 20. 
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Fig. 9. — The gravitational radiation waveforms for an observer situated along the rotation axis (z 
axis) in geometrical units as a function of time. (a) The + polarization. (b) The x polarization 
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Fig. 10. — The energy emitted in gravitational waves for the P and N runs. The calculation was 
started at t = —25 when the runs had about the sanie relative separation. 
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Fig. 11. — Coordinate rest mass density contours for the runs. Time increases to the right. The 
length scale is in R^,. The contours are taken at the initial time step, t = —40, t = —25 (beginning 
of dynamical instabihty) and t = — 15 (a = 2i?*). The rest mass density is in units of M/R^ and 
the contours are logarithmic with a spacing of 2.3 starting at 1 
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Fig. 12. — Coordinate rest mass density contours for the runs. Time increases to the right. The 
length scale is in i?*. The contours are taken at t = — 7 (peak in Fig 7c), t = O (minimal separation), 
t = W and t = 20. The rest mass density is in units of M/Rf and the contours are logarithmic 
with a spacing of 2.3 starting at 1 
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Fig. 13. — Coordinate rest mass density (p^,) in units of (M/R^) for the final configuration [t = 20) 
on the X — z plane. The contours are logarithmic with a spacing of 2.1 starting at 1. The length 
scale is in R^ 



